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Abstract 


A theoretical model of a laminar diffusion flame at the 
leading edge of a fuel plate in a forced convective flow is pre- 
sented and solved numerically to study the flame stabilization and 
blowoff phenomena. The system of governing equations yore of 
the two-dimensional Navier-Stokes’ momentum, energy and species 
equations with a one-step overall chemical reaction and second- 
order, finite rate Arrhenius kinetics. The computation is per- 
formed over a wide range of Damkohler numbers. For large Damkohler 
numbers, envelope flames are found to exist where the computed fuel 
evaporation rate, the flame stand-off distance and the velocity 
profiles show certain similitude. As Damkohler number is lowered, a 
transition to open-tip flame takes place where the flame becomes 
stabilized on the sides of the fuel plate. Further decreasing of 
the Damkohler number pushes the diffusion flame downstream out of 
the leading edge region. In this paper, the flame structures of 
the envelope and the open-tip flames are presented together with a 
description of the transition sequence. The implication of this 


work to downstream boundary layer combustion is also discussed. 
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CHAPTER I 


INTRODUCTION 


1.1. The Need for an Elliptic System of Equations 

The burning of a thin fuel plate placed in a parallel oxidi- 
zing stream offers a farthcok otiar ls to illustrate the interaction 
between aerodynamics and heterogeneous combustion [1, 2]. In the 
limit when chemical kinetic rates were much faster than those for 
heat and mass transfer, Emmons [3], using boundary layer theory, 
showed that the evaporation (pyrolysis) rate of the condensed fuel 
could be expressed as a function of the Spalding mass transfer 
number alone. Because of the simplicity of the result, this type 
of analysis has subsequently been adopted in many applied problems 
[4]. 

The use of infinite fast kinetics pre-assumes the existence 
of the diffusion flame. When the question of flame blowoff or 
extinction is at issue, the analysis has to include finite-rate 
chemical kinetics at the leading flame front region. This invali- 
dates the boundary layer approximation because of the reasons to be 
discussed. 

There are two main reasons why boundary layer approximation 


is inadequate in the leading flame front region: one consideration 


has to do with combustion, the other has to do with fluid 
mechanics. The former consideration will be emphasized in this 
discussion. For the problem concerned, we assume that flow 
separation and recirculation do not occur, the temperature of 
incoming oxidizer flow is low and the activation energy for 
chemical reaction is high. Under the circumstances, the initiation 
of chemical reaction in the leading flame front region requires an 
efficient heat conduction from the downstream hot product zone to 
the cold upstream. In other words, without heat conduction, the 
flame will be unable to be stabilized. Therefore, for analyzing 
flame stabilization and blowoff, streamwise heat conduction has to 
be included. This renders the partial differential equations el- 
liptic in nature. 

We also ane gbss at this problem from a different angle. The 
mathematical problems of flame blowoff and extinction are charac- 
terized by the existence of multiple steady-state solutions [5]. 
If streamwise diffusions (momentum, heat and mass) are neglected in 
the governing equations, the system is not capable of producing 
multiple solutions. A unique solution is determined by the up- 
stream conditions and the boundary conditions at the wall and the 
free stream, because the governing boundary layer equations are 
parabolic in type and can be solved by the space-marching technique 
[6]. One exception is the stagnation-point boundary layer which 


yields an ordinary differential equation after transformation and 


in this special case, streamwise diffusion is retained. 

In addition to the above combustion consideration, the accu- 
racy from fluid mechanical analysis also requires the inclusion of 
streamwise diffusion terms. The characteristic length in the 
leading flame front zone is given by the thermal diffusion length 
which normally is very small. The streamwise gradients, therefore, 
can be comparable to the cross-stream ones. Additionally, for the 
particular problem considered in this paper, the leading edge of 
the flat plate has been well known for its necessity of elliptic 


treatment [7-9]. 


1.2. Related Literature 

Probably because of the mathematical difficulty associated 
with the solution of elliptic type partial differential equations 
there are not too many theoretical treatment on stabilization and 
blowoff for multi-dimensional diffusion flames. The works cited 
below are all solved by computational techniques. 

To the best knowledge of the authors, the first successful 
analysis addressing the extinction aspect of a two-dimensional 
diffusion flame is the work of Frey and T’ien [10]. Their work 
offered a diffusion flame spread model over a thin fuel in an 
opposed flow. The emphasis is concentrated on the energy and 
ereties equations using a simplified assumption of constant 


density, constant properties and uniform velocity. To include the 


fluid mechanics-combustion interactive effect, reference 11 gives 
an exact treatment of Navier-Stokes’ momentum equations, but is for 
a premixed flame stabilized in a tube. References 12 and 13 treat 
a diffusion flame stabilized over the surface of a solid fuel. 
These works include complete momentum equations, but are done 
dimensionally. Reference 12 presents flame structures for both 
forced and buoyancy-induced flows. Reference 13 studies flame 
front positions as a function of convective velocity. Due to em- 
ploying an inert plate in front of the pyrolysis zone, the model is 
not applicable in the vicinity of the leading edge of the fuel 


plate. 


1.3. Problem Definition 

The specific problem to be treated in this work can be illus- 
trated by Fig. 1. A two-dimensional, thin and flat fuel plate is 
placed in a forced oxidizer flow which has uniform velocity and 
temperature profiles upstream. Under proper conditions, an envelope 
diffusion flame is established ahead of the leading edge of the 
fuel plate as indicated in Fig. 1. The flame transfers heat in both 
x (streamwise) and y (cross-stream) directions to preheat the 
incoming cold oxidizer and to vaporize the fuel. Similarly, the 
Mass transfers in both x and y directions supply the balanced 
amount of oxidizer and fuel vapors to the flame reaction zone. The 


purpose of this theoretical work is to determine under what 


conditions an envelope flame can exist, its structure and to inves- 
tigate other flame configurations when the envelope flame is no 
longer possible, such as blowoff. 

Because of symmetry with respect to the x-axis, the domain to 
be studied will be reduced to the half plane o ¢ y<e. The coor- 
dinate system is such that the origin (x=y=0) is always located at 
the leading edge of the fuel plate. No burnout of the fuel sample 


is included. 


CHAPTER II 


MATHEMATICAL MODEL 


2.1. Introduction 

Since the problem has been defined, the combustion model for 
the flame stabilization zone can now be set up. It consists of 
conservation equations for mass, momentum, energy and species in 
differential form. An equation of state, the variation of vis- 
cosity with temperature and a number of other assumptions are used 
as the linkages between properties. An assumed chemical reaction 
rate expression is also included in this system. Finally, a set of 
known boundary conditions is specified to complete the mathematical 
description for this problem. 

Before the actual computation, a normalization procedure is 
performed in the model by using proper length scale and reference 
quantities. The purpose is to seek out the nondimensional para- 
meters, which appear directly as the coefficients of unit order 
terms in equations, for conditions of interest and to evaluate the 
relative importance of various terms. In addition to simplifying 
the equations, the nondimensional procedure also leads to simi- 
larity considerations as well as to reducing the cost and the 


effort for numerical computations. 


2.2. Conservation Equations in Dimensional Form 

The conservation equations for multicomponent reacting mix- 
ture are well established in many references [1, 2]. Before 
writing them down, the following assumptions are made: 


(a) the flow is steady, two dimensional and laminar 
with Stokes’ hypothesis, 


(bo) radiation heat transfer and the influence of 
body forces are negligible, 


(c). viscous dissipation and compressive work are 


neglected owing to the low speed combustion 
problem, 


(d) the ideal gas law is applicable to the gas 
mixture with constant and equal specific heats, 
equal binary diffusion coefficients, constant 
Prandtl and Lewis numbers and constant value of 
PH. 

(e) the average molecular weight is constant. 


Based on the above assumptions, the mixture equations, with 


the coordinate system defined in Fig. 1, are: 


continuity equation: 


ou) | 3(ev) _ 4 (2.1) 
ar Te aa g 
ax dy 


momentum equations: 


ce act ices ue Sa in set 
Gly 2 Rep S Say dl Tabee ss cay MR) ce pine eee 
ay = pa eel pr (22) 
= ele an tae aly 


ser .5PF. 2s -s By (a8 , 37 
GEOR sntoow TewaquPaiinc toy uetennisy setae 
dues skdv 2,0u OV 
+ gylulesy - 35m + gy)! (2.3) 
energy equation: 
—- oT <-_ aT a, c3T of aT eee 
ee — ee — ee — er —— 25 = Seer} 2.4 
PucEsg + PYCEaS 5x XS) + ay nal + Q (2.4) 
fuel species equation: 
eae ony: ~-9Y. _80¥ rer —< 
eB, 


FE AWK AEWA ER ie 
D—=) + ay ODF) + Wp 


Oxidizer species equations: 


aY oY aY oY = 
ae LD, aes oO Fb ig = eG bic am eo ° 
Ou ane + Sieg = 5x PD) + Sten? + Wo (2.6) 


equation of state: 


p = ot, Cp 


and the viscosity variation with temperature is taken to be 


(2.8) 


rl 
iT} 

a 
‘“d 


Since the convective flow velocity in this problem is much 
slower than the speed of sound, the ratio of the pressure 
difference to the absolute pressure throughout the region of 
interest is much Nererchad one which implies that the dependence of 
density on pressure can be ignored. Therefore, the equation (2.7) 


can be rewritten as 


aoe (2.9) 


2.3. Chemical Kinetics 
In the present study, a one-step overall chemical reaction is 


proposed. Gaseous fuel reacts with oxidizer to form product and 


10 


release heat, or 
1 Kg [Fuel] + f Kg [Oxidizer] 
—-= (1+£) Kg [Product] + Q KJ [Heat] (2.10) 


where f and Q represent stoichiometric oxidizer/fuel mass ratio and 
heat of combustion per unit mass of fuel respectively. The re- 
lationship between mass and energy source terms is 


bp = up/f = -Q/Q Pepey 


The rate equation for fuel is assumed to be concentration 
dependency of second order and temperature dependency from 


Arrhenius’ Law, i.e., 


Op = '-Bp2¥pYq EXP (-E/R°T), gisey 


where B and E are the frequency factor and the activation energy 


respectively. 


2.4. Boundary Conditions in Dimensional Form 

To complete the model, it is necessary to specify the 
boundary conditions. From the configuration of Fig. 1, due to the 
symmetry with respect to y=o line, the computational domain is 


given by 


ats | 


where Xnin 224 Ymax #7e chosen to be distance large enough so that 
the influence of the fuel plate and the flame are negligible. The 


value of x (downstream) is chosen where the streamwise gradi- 


max 
ents are expected to become small enough. 


The boundary conditions are: 


at x = x13, (upstream): 


SPROUL, PFS OF Whe ree lyse cto, pyepisnye (2.13) 
AS es ee (downstream): 

jeeovesT SE tO | (2.14) 
eee bh? [er ery ee ’ ; 


aitpyt=cyiet (far from the stagnation stream line and the fuel 


plate): 
au =U Pe, TetndThy Wie aic0)1 Yeiter ¥ (2.15) 
a eae cris Ba ay aac? ae POS aq SOF: 


where the zero derivative of v results from the application of 
continuity equation over there. Along y = 0, where the plate 


leading edge, x = 0, separates the gas phase stagnation streamline 


12 


from the solid fuel plate, we have 


for Xin< x < o (on the stagnation streamline) 


=- = oY oY 

du = oT F 0 

—w = Vis —s = —» =- 0, —s = 0 ’ 
ay oF 0, ay 0, ay EE ’ (2.16) 


due to symmetry, and 


for 0 < 2 < X,4, (on the plate): 


u = 0; T S Tye 
= Stee 
mY oy =m + eD—sl.,, 
= nit 
my ow ? PDs ly, 
Bre op a it 22 59 
= iF lw af alle 
= Pere py 
= 0 f§ =< 0 
z mech e Neaas (2.17) - 


where m = pyv, is the unknown mass burning rate on the solid 


surface. Since the solid fuel surface temperature is maintained 


13 


constant and greater than the ambient value, the last condition in 
Eq. (2.17) is imposed to prevent negative burning rate. This 
artifice does not change the essential feature of the problem and 


can be replaced by using a surface pyrolysis law. 


PA Ye Dimensionless Governing Equations and Boundary Conditions 
A length scale and the reference state should be specified 

properly prior to the normalization procedure. The characteristic 

length chosen for nondimensionalization is the thermal diffusion 


length given by 


‘nae (2.18) 


where a is thermal diffusivity. U, is free stream velocity and 
superscript * denotes quantities evaluated at the reference tem- 
perature je which, in this study, is chosen to be the arithemetic 
average of the ambient and the adiabatic flame temperatures. The 
correct choice of the thermal diffusion length as the reference 
scale will be reflected by the simplicity of the dimensionless 
equations, the magnitude of the computational domain and the 


similarity of the results. 


The decision to evaluate the reference thermal properties at Ta 


is based on the consideration that many of these quantities vary by 


14 


one order of magnitude within the flame due to the temperature 
difference, so evaluation at the average temperature is a better 
representation of the magnitude. However, the temperature itself 
is normalized by T,, rather than T°. This will be discussed later 
in Section 4.1. (Parameters). 


The nondimensional variables are defined as:- 


Pee ae aero i fee 
Oh o*! in Mh o*y2 ’ O*u ? itis? 
k* =  (6/Ue) 
Da = (~73)8 = (T/p*E)’ 
Cut 
p-o2 
me tiesto Pas rae 
T = 7" Q EG ae E=Rqy LG T° epee ea 
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Based on the length scale and reference quantities just 
defined, a normalization procedure is carried out through the 
model. The dimensionless governing equations and boundary 


conditions are presented as follows: 


continuity equation: 


“1 Bh (2520) 


gu gu, 9P , 9,u ,5,au _ 2,3u , av 
Rome wuacx ax Re\iax , 3'9x “ay .¥ 
Ta Yuu du *8 ov Q421) 
* Ene ome Vee 
av Gu ersD (hee ey 80 ey 
eee Sal taryee MaTyie oC Une: jyrg: 


u,.~oVv eeu ov 
: nae sable tonu eg (222) 
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where Re is the Reynolds number defined by 


energy equation: 


oT 


fuel species equation: 


oY oY oY oY. 


1 
+ OV 


: me 2 a 8 
“he Jy RePrLe 9x4 5x) Y use SH 


ox 


+ Way 


oxidizer species equation: 


oY 3 
O Yo oY 
pu Oy ay ee ae SCA geen 
Pa! dy “ Reprie ax) max’ * syumoede 
+ fw 


F’ 


(2.29) 


(2.24) 


PAS 


(2726) 


Li 


where 


= -Dap“YW¥, EXP (-E/T) (2.27) 


Or F-0 


is the nondimensional fuel reaction rate, 


equation of state: 


Ql ee (2.28) 


and the viscosity variation with temperature is taken to be 


nei ae (2.29) 


Using Eq. (2.18), the Peclet number is found to be unity, i.e., 
RePr=1. Therefore, Eqs. (2.24)-(2.26) are greatly simplified. 
Since Re=1/Pr and Prandtl number for gases is of the order unity, 
the flow problem near the flame front is of low Reynolds number 


type. 


18 


The boundary conditions are 


at x= Snin’ 


F ei ty en’ 
at x ™ Inex 
ul ravi STE oa Ome 
X x ax ox gx 
at y = Ynax 
Ba ilcete ucbtael «| uae Oy 684 NO 


(2.30) 


(2.31) 


(2.32) 
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oY 
a F 
By i ey Cae tage hee ee 
< * < “Tmax 
’ m = Pag! T = Ty? 
; Thad bd 
= M+ RaprLe”aylw 
' 1 he 
~ RePrLe™ ay lw 
1_,,22) .¢ oT 
RePr'oy W * dy 
thet T 
+f oY 


(2.33) 


(2.34) 


CHAPTER III 


NUMERICAL SCHEME 


3.1. Introduction 

The present model is solved numerically. The method of 
deriving the discretization equations for given differential 
€quations is to use the idea of control-volume, The most 
attractive feature of the control-volume formulation is that the 
resulting solution implies that the integral conservation of 
quantities, such as mass, momentum, energy and species is exactly 
Satisfied over any group of control volume. So it directly leads 
itself to the physical interpretation. 

The arene? primitive variables, rather than the stream 
function/vorticity formulation, has been chosen in this work. One 
of the difficulties in using the primitive variables in solving the 
Navier-Stokes’ equations is the lack of an explicit equation for 
pressure. A special treatment originated from the SIMPLE procedure 
developed by Patankar and Spalding [14] is introduced to resolve 
the pressure gradients in the momentum equations. Staggered grid 
system is utilized in the numerical scheme. Since the differential 
equations are of elliptic type and coupled together, an iterative 


procedure is needed in the computation. 
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pia 


In this chapter, the algorithm for the numerical scheme 
including the description of grid system, the finite difference 
equations, the specified pressure correction equation, the 
computational domain, and the boundary conditions is developed 
first. Next, a sequential procedure to apply the alogrithm is 


presented. 


3.2. Grid System 

The calculation domain is divided into a number of 
nonoverlapping control volumes such that there is one control 
volume surrounding each grid point. This pattern of two dimen- 
sional grids is shown in Fig. 2. The variables, except the veloc- 
ity components u and v are stored in main grid points. While nu and 
v, represented by short arrows, are placed on half-way between the 
two adjacent grid points in Fer) y directions respectively. 

The grid system now can be categorized into main cell, u cell 
and v cell according to the control volumes in relation to each 
assigned location of variables. They are shown in Fig. 3, 4, and 5 
respectively. The subscript P in these figures denotes the 
concerned point, whereas E, W, N, and S mark the neighbor points 
surrounding P. In main cell, it is easy to see that the 
discretized continuity equation is applied by the differences of 
adjacent velocity components, and the convective terms for 
variables are defined just at the faces of control-volume. As to 


the streamwise velocity u, the control volume, defined as u cell, 
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is staggered half-way in x direction with respect to main cell. The 
choice of v cell for normal velocity v has similar behavior as u 


cell except that the staggering movement is along y direction. 


3.3. Finite Difference Equations 

The finite difference equation for each variable is obtained 
by integrating the differential equation over the relative 
computational cell, associating with specified interpolation 
applied on the interface between the variables of two adjacent grid 
points. The detailed derivation of finite difference equations is 
carried out in Appendix A. The resulting finite difference 


equations are: 


continuity equation: 


CU t=) fin ue Coe me eee ee, 


momentum equations: 


a ein PE u u 

apup = va nb * Sa + Aaw'Pp - Ph), (3.2) 
V V V 

Sp’ putetanpYnpet Se chwApeGpperiep, ). Ooty 
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The energy, fuel and oxidizer species equations can be expressed as 


following general form of 


= o 
p°p = cab nb or Sas (3.4) 


where represents temperature, fuel and oxidizer mass fractions. 
The summations are applied over E, W, N, and S points, the values 
of dependent variables u, v, >, and p are evaluated at (n + 1)th 
iteration and the expressions of coefficients such as A’s, a’s, S’s 
and C’s are given in Appendix A. 

The appearance of Equ. (3.1)-(3.4) seems to be linear, but it 
is not because the coefficients are also the functions of dependent 


variable itself. 


3.4. Pressure Correction Equation 

The momentum equations can be solved only when the pressure 
field is given or somehow estimated. An equation for pressure is 
needed to resolve the problem. The thinking of SIMPLE procedure is 
that unless the correct pressure field is employed, the resulting 
velocity will not satisfy the continuity equation. Therefore, 2a 
pressure correction equation is derived by linking the continuity 
Cquation. The derivation is outlined next. 


Since the pressure field is unknown in the beginning, a p* is 
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guessed in momentum equations. An imperfect velocity field based 
on guessed p* will be denoted by au’ and v*. This crossed velocity 


field is resulted from the following equations. 


us+ u u + + u 
V+ V o+ Vv + +,,V 


The continuity equation generally will not be satisfied by 


+ + 


and v into it, instead a net non-zero mass source mn 


employing u p 


will be generated where 


+ + 
P Ee WE NP SS¢ Se 


Since the mass source is originated from pr’, the next step is 
to find a way of improving the guessed p* such that the resulting 
crossed velocity will progressively get closer to satisfy the 


Continuity equation. First, a correct pressure is proposed that 


(3.8) 
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where Dp is called the pressure correction. Suppose that the true 
velocity components u and v respond to this pressure change in the 


following ways: 


re + ' ’ u u 
Up = Up + (Py - Pp) (Ag, /ap)s , (3.9) 


~ a yes ' Vv Vv 
Ves a Yop te (Dp Py) (A, 5/ap)- (3.10) 


Then, substituteall the velocitycomponentsinto the continuity 
equation using the above velocity-correction formulas (Eqs. (3.9) 
and (3.10)), we can obtain, after rearrangement, the following dif- 


ference equation for p: 


Pius Pit Pus Pit Pp 
EPE a anew + anPy + asPo + Sa p92 be | 
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where 


p u u 
w= CwlAaw/@plp. 


fu 
0 


Dp ei aatd 
CslAns/apl ps 


Po ep Pp Pp p 
ap = as + any + an + ag. 
Sapa 
Sa = Mp - 


Once p is obtained from Eq. (3.11), the velocities and 


pressure will be updated by p’ through Eqs. (3.8)-(3.10). 

In the computation, as soon as the numerical results meet the 
criteria of convergence, which will be discussed in Appendix B, the 
value of aoe) will come out to be practically zero for all control 
volumes. Therefore, py’ =o at all grid points will be acceptable 


solutions of Eq. (3.11) and the crossed velocities and pressure 


will be the correct velocities and pressure. 


3.5. Computational Domain and Grid Size Distribution 


Except for the fuel plate and the stagnation streamline, 
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which are located at y = 0, the other three boundaries, Sern ana y 
and Yiax @re in the gas phase where the appropriate locations have 
to be determined. Mathematically, these values should either be -< 
(xnin) OF +™ (xya5- Vmax)? DUt for numerical analysis, appropriate 
finite values are needed for efficient computation. 

The absolute magnitudes of x;, and y,,, are consideredlarge 
enough when the boundary conditions, Eqs. (2.30) and (2.32), are 
satisfied within a given margin and the corresponding gradients 
become vanishingly small at the same location. One may intuitively 
expect that this occurs at a distance several times of the flame 
radius which, in itself, is several times of the reference thermal 
length. In practice, a numerical experiment is required for 
optimization. In this study, z,;, = -20 and Ymax = 60 have been 
found satisfactory. 

The location aax for the downstream boundary needs somewhat 
different consideration. Far downstream from both the leading 
edges of the flame and the fuel plate, it is expected that boundary 
layer approximation will prevail such that d/dy >> d/dx. Since the 
derivatives with respect to y are expected to be of the order of 
unity, the boundary conditions, Eq. (2.31) for 0/dx = o are adopted 


@s an approximation applied at a location (x ) where 


max 


boundary layer approximation is valid. The question is then what 


is the minimum acceptable x for both numerical accuracy and 


max 


computational economy. 
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A series of numerical experiments were performed for nat 


equal to 15, 30 and 40, and comparison of results in the domain of 
interest (x < 10) was made. The difference between the results of 


Jwexct 30 and «x éz0" 40 is negligible. The comparison between 


x = 15 and x = 40 is shown in Fig. 6. While their dif- 


maz max 


ferences are not large, they are finite, due mostly on the dis- 
placement of the flame position. However, we have decided to use 
Imax * 15 and to Droste t this margin of inaccurarcy in the sub- 
sequent parametric study for the sake of computational economy. 
More discussion of the evolution to boundary layer can be found in 
Section 4.5. 


A variable grid size distribution is adopted in this study. 


The grid increments are designed according to the following 


formula: 

u u m 
(Ax)mog mAl Aragon +4 1)", 
(AyJivews af Ay) Fea (Gane 


where (Ax),” and (Ay), % are the initial grid lengths surrounding 


the plate leading edge and they are also the smallest sizes in x 
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and y directions. Grid sizes increase from the origin toward the 
upstream (-x), the downstream (+x) and the outside free stream 
(+y). The rate of expansion near the origin is initially slow, but 
accelerates quickly away from the origin. 

In this work, (Ax)," and (Ay),", are taken to be 0.2. Using 
this arrangement of grid spacing, and the computation domain just 
discussed, the mesh points are 49 x 38. The grid distribution is 
presented in Fig. 7. An accuracy check on the grid size and 


its distribution will be presented in Appendix B. 


3.6. Boundary Conditions 

The boundaries of the computational domain for the flow field 
are placed on the mid-way between the main grids so that the 
velocities normal to the control-volume faces (coincident with the 
boundaries) are immediately prescribed. The linear interpolation is 
applied on the boundaries for the other variables located on both 
sides of each boundary. 

In the SIMPLE alogrithm, it solves the pressure correction 


) rather than the absolute pressure p itself, therefore, the 


P 
boundary conditions for p’ need to be specified. 


The arrangements of boundary conditions are as follows: 
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(a) at x = X55: 


The boundary conditions are 


Uti joan 
V'Cie on ea 
D1, J) 18 1) 
Yo(1,J) z 0 


Since the velocity u on the boundary is given which implies 
there is no correction for itself over there. To apply Eq. (3.9), 
ur is equal tou that leads the normal gradient of Pp to be zero. 


Therefore, the boundary condition of p is 


piott;J)qestp? cae t 


(b) stex =x 


A special treatment of u to specify 
p on downstream boundary is adopted from 
Ref. 11. The value of u* computed in the 
last iteration, by using the condition 
that the normal gradient of mass flux at 
this boundary vanishes, is taken as the 


boundary value for u at present iteration. 


31 


As a result, it imposes zero normal 
gradient condition on p’ . The boundary 


conditions are: 


u(MI,J) = u*(MI-1,J) 
v(MI+1,J) = v(MI,J) 


T(MI+1,J7) 


T(MI,J) 


Y¥5(MI+1,J) = Yo(MI,J) 


YQ (MI+1,J) YO (MI,J) 


p'(MI+1,J) BeeMIs oh) 


i4= MI-1 i4s MI 


where the superscript + denotes the value in the last iteration. 


cimoet y= yo: 


Since the normal gradient of v is 
zero, the same procedure as used for the 


boundary at x = x1,, is applied for v and 


p’. 


The boundary conditions are 


GOLeMo el yaad 
ViryMd) = vel M1) 
TOGIM IST) = 1 


Yo(I,MJ+1) 0 


Yo 


DG GM) 


Yo(I,MJ+1) 


Otol Mo + 1h) 


32 


(d) at y = 0: 
The velocity v normal to y = o is specified either by the 


symmetrical condition with respect to the stagnation streamline or 
by the coupled relation with the condition as gas-phase, therefore, 


it results that the zero normal gradient is applied forp. The 


boundary conditions are: 


(dol)ox3n |< 3 < 0: 


UO CI4r) em ae 2 

V (Tp t-) a9 

TED ps) eel el, 2) 
Yp(I,1) = ¥,(1I,2) 
Yo(t,2) 
Piatt 2) 


Yo Te) 


pitty) 


(di 2) agome XS XInax: 


u(I,1) = -u(I,2) 
T(I,1) = 2T. a ies hi EA 


m 


w p.jv(I,1) 


Be GL se) = op] ae 


na ° 
WwW Uw ag Uw _ ma 
ra Repruedy pit!) ae ve ca (RePrLedy 2)¥p(1,<) 


29) 


m : 
W Uw uw m 
(aya RGESTEay mS aceyleht 2 (RePrLeby iy PET ee 


siete be 
Since the fuel plate involves the leading edge, the 


expression of velocity at x = y = 0 needs a special treatment. When 
the discontinuity is of the functional difference form, lim v(o ) = 


lim v(o*) it can be evaluated by using [15] 


Mi geteny iwi! 2 (MH, 2) = TM ND 


IE ae eae ie i oe 


where (MH,1) is the location of leading edge. 


3.7. Computational Sequence: 


The order of computational execution for a particular 


iteration is: 


1. Specify the values of all the variables including 
p for each grid point in the flow field. These 
data either come from the results of previous 
iteration or initial guesses. 
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2. Solve the momentum equations, Eqs. (3.5) and (3.6) 
to obtain u* and v*. 


3. Solve the p’ equation. 
4. Calculate p from Eq. (3.8) by adding p’ to p. 


5. Update u and v from their crossed values by using 
Eqs. (3.9) and (3.10). 


6. Solve the discretization equation for other 
such as temperature and species. 


7. Treat the corrected pressure p as new guessed p* 
and update all the variables, return to step 2 and 
repeat the whole procedure until a converged 
solution is obtained. 

During the above sequence of computation, the dependent var- 


iables and p’ are solved via their corresponding finite difference 


@€quations which have an identical form given by 


The formula itself is nonlinear, therefore, an iterative method is 
tsed to solve the problem. The coefficients in the equation are 
first linearized using the values evaluated from the previous 
iteration. Next, the linearized algebraic equation follows a line- 
by-line iteration which usually accelerates the rate of con- 
vergence. The line-by-line iteration is commonly known as the ADI 
(alternating direction implicit) method which is first developed by 


Peaceman and Rachford [16]. The complete iteration consists of a 
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first half in the x direction followed by a second half inthe y 
direction through successive employment of the tridiagonal-matrix 


(TDMA) at each sweep. The iteration proceeds in the x direction 


that is 
1 1 | 
aie Pict n+s n+> n n 
apo, ~ = apdn ~ + ayo, (any + 45% 
+ SQ) Pee (212) 


where the superscript n + 1/2 marks the solution which is going to 
be resulted from this sweep, whereas o's and ce are known values 
from the last (nth) iteration or from the initial guesses. The 


followed iteration in y direction is 


Neto. n+1 n+1 n+s n+s 
Qap?p = ayoy Ot AOS + (apop * + ayo, ey en 
+ Sq) 


where superscript n + 1 denotes the unknown values. The appearance 
of Eqs. (3.12) and (3.13) is of tridiagonal form, the TDMA is used 
for solving each line iteration. 

Since the algorithm and the method of solution have been 
introduced, the final concern for the numerical scheme is conver- 
gence. To guarantee convergence, the algebraic equations must 


satisfy certain conditions from the theory of linear algebra [17]. 
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The difference equations can be written in following general form 


when the equation is linear, (a’s and b’s are constant), the 


formula will give convergent results if: 


Sou eae 1 (3.14) 
s%iay lis! < | 


for each i, and 


jeiedtsiglion & 


for at least one i. 


In the present derivation of finite difference equations at 
each iteration the coefficients are treated as constant and the 
Coefficients indeed meet the conditions given by Eqs. (3.13) and 
(3.14). Therefore, the convergence is ensured. However, the co- 
efficients in the finite difference equations change from iteration 
to iteration, sometimes divergence may occur. To avoid this situa- 
tion, an underrelaxation factor a. is applied which can be ex- 


Pp 


pressed as 


oH 


where $* are used to update the coefficients in the next 
iteration. 
The typical considerations of convergent criterion and the 
convergence time will be discussed and presented in Appendix B. 
The results to be presented in the next chapter are obtained 
using a computer program written by the author based on the SIMPLE 


procedure just mentioned. 


CHAPTER IV 
RESULTS AND DISCUSSIONS 


4.1. Parameters 

There are ten nondimensional parameters, as shown in Table 
4.1., derived from the normalization procedure for the system. We 
Mentioned previously that the thermodynamic and transport pro- 
perties are normalized by their corresponding quantities evaluated 
at the reference temperature T* with the exception of the tem- 
perature itself which is nondimensionalized by the ambient tempera- 
ture T.. If we choose the reference temperature for property 
evaluation to be Ty, instead of T°, the number of nondimensional 
Parameters will be nine rather than ten. However, this shortens 
the thermal length (Eq. 2.18) and consequently the computational 
domain in terms of nondimensional distance has to increase accor- 
ding to the ratio of To sre By choosing the average temperature as 
the reference temperature, this hidden parameter, Tare is brought 
to the surface and the reference length scale becomes a more rea- 
listic one. 

Since the nondimensional parameters are made of groups of 
physical properties, they are needed to be specified initially. The 
values of properties are listed in Table 4.2. It should be men- 


tioned that the data used in gas-phase chemical reaction such as 
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TABLE. .4.1 


Nondimensional Parameters 


Symbol Parameter Group Value 
Pr v/a | 1 
Le . a/D 1 
Fie ee PANT 
7 Neate Oh 4,33 
8 ape A. a) caPae pi 
f 3 an 
Q Eisele tt 116.8 
L Ry cst 4.33 
E E/R°T) 54.20 


Da (k*/c,U )B Vanlas. 


Name 


Ambient Temperature 
Reference Temperature 
Density (reference) 
Viscosity (reference) 


Thermal Diffusivity 
(reference) 


Specific Heat 

Frequency Factor 
Activation Energy 
Heat of Reaction 


Evaporation Temperatu 
of Fuel 


Latent Heat 


Free Stream Velocity 
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TABLE 4.2 


Property Values 


1g (= T 


Value 


300 
hen OD 
0.2/1 


4.93x10 5 


2°. 58 x0r 


Meus ial 


1.58x10°2 


35x07 


47,400 


650 


I Smee) 


variable 


5 


4 


Unit 


KJ/Kg°K 
m?/Kg"S 
KJ/Kg-mole 


KJ/Kg-fuel 
K 


KI/Kg 


m/S 
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average specific heat, adiabatic flame temperature, stoichiometric 
mass ratio, heat of combustion, pre-exponential factor and 
activation energy are adopted from Lee’s work [11] where the fuel 
is taken as methane (CH,). 

The effect of varying Damkohler number is investigated in 
this study. The Damkohler number, Da, defined by Eq. (2.19) is 
proportional to the ratio of gas resident time in one thermal 
diffusion length (5/U,) to chemical reaction time (1/p"B). 
Finite-rate chemical kinetic effect is, therefore, expected as 
Damkohler number becomes small enough. For a given fuel and orxi- 
dizer, small Damkohler number can be achieved by increasing the 
forced free stream velocity U,. 

Using the values of nondimensional parameters listed in Table 
4.1., the computation is performed for a wide range of Da from 
5.32 x 10° to 4.50 x 107 which is categorized into 7 cases as shown 
in Table 4.3. For Da>7.03 x 105 (Cases 1-3), envelope diffusion 
flames surrounding the leading edge of the fuel plate are found. 
For 5.82 x 10° < Da < 7.03 x 109, a transition between envelope and 
Open-tip flames occurs (Cases 4 and 5). For 5.32 x 10° cea 15.82 
x 10° (Cases 6 and 7), open-tip flames, where the flame front is 
retreated from the upstream to the sides of the fuel plate, exist. 
For Da < 5.32 x 105, the diffusion flame is blown off the computa- 


tional domain of the leading-edge region (x < 10). 


42 


TABLE 4.3 


Parametric Studies 


Case Damkohler Number (Da) Type of Flame 
1 4.05x107 E* 
2 2.81x10° E 
3 7.03x10> E 
4 6.38x10> T 
5 6.08x10> T 
6 Baise Tae 0 
7 5.32x10> O 


* Es: Envelope Flame 
T: Transitional Flame 
O: Open-tip Flame 
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Computed flame front structures for the envelope, tran- 
sitional and open-tip flames will first be presented in the next 
three sections. A consideration of flame similitude for envelope 
flames and a discussion of the implication of leading-edge flames 


to boundary layer combustion will then follow. 


4.2. Envelope Diffusion Flames 

The numerical computations show that for Damkohler number 
between 7.03 x 10° and 4.50 x 107, envelope flames exist 
surrounding the leading edge of the fuel plate. Although no 
computation has been made for Da greater than 4.50 x 107, it is 
expected that envelope flames should persist into the larger 
Damkohler numbers for the acchentt esi system we are dealing with. 

Comparison of the computed results for the three envelope 
flames, Cases 1-3, indicates that most of their essential features 
are similar to each other. Therefore, the narration for the 
detailed flame structure will be done only for Da = 2.81 x 10° 
(Case 2) using Figs. 8 to 13. Corresponding flame profiles for 
Case 1 and Case 3 are presented id Figs. 23-32, for those readers 
who need more details. Additional discussion on the similarity and 
the difference between the three cases will be made later, in 
Section 4.5. 

Fig. 8 is the distribution of isotherms in the region of 


interest. An envelope flame is situated upstream of the leading 
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edge of the fuel plate. Along the stagnation streamline (y =0), 
the maximum temperature is located approximately two thermal 
lengths from the leading edge. In the vertical direction the flame 
stand-off distance is greater because of the blowing effect of fuel 


vapor at the wall. For example, T along x = 0 occurs at y = 4. 


max 
The magnitudes of the flame stand-off distance and the flame 
thickness suggest that the choice of the reference length is 
reasonable. The radius of flame curvature is comparable in magni- 
tude to the flame thickness, confirming the two-dimensional feature 
of the flame front region. Note that the highest temperature in 
the field (T > 6.9) is located off the center line y = 0, its 
precise position is sensitive to the values of the Damkohler 
Dumbers, as shown in Fig. 15. Fig. 9. gives an alternative presen- 
tation of the temperature profiles. Along the symmetry line y = 0, 
Lo occurs at x = -2.40, its magnitude is only slightly lower 
than the highest temperature in the field. The temperature drop 
Near the stagnation line between x = -2.40 and x =o is due to the 
quenching effect of the cold fuel plate. This figure also indi- 
Cates that the temperature profiles are much fuller on the fuel 
side than those on the oxidizer side. 

The fuel and the oxidizer mass fraction distributions are 
illustrated in Fig. 10. Since the flame encloses the leading edge 


of the fuel plate, only a small amount of fuel vapor is able to 


penetrate across the flame because of high reactivity there. Simi- 
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lar situation is applied for the oxidizer. Approximately, we can 
speak of a fuel and a oxidizer sides. Near the flame tip, the 
large fuel concentration difference produces a strong mass dif- 
fusion opposed to the convective flow which establishes the reac- 
tion zone ahead of the fuel plate. 

Fig. 11 is the distribution of nondimensional reaction rate 
of fuel, OF = Dap*YpY,EXP(-E/T). The most vigorous reaction occurs 
around the apex of the flame tip, but is slightly off the 
stagnation streamline. Unlike the predication by the thin-flame- 
sheet theory, the highest reactivity is not coinciding with the 
highest temperature region (T > 6.9). This behavior appears to be 
Originated from the effect of finite-rate chemical kinetics and 
from the definition of WE the reaction rate not only depends on 
temperature but also on the product Yr Y,. Similar observations 
have been found for the transitional (Figs. 35 and 40) and the 
Open-tip (Fig. 18) flames, as well as in Refs. 12 and 18. The 
strong reactivity zone around the flame tip serves as the stabi- 
lizer or the igniter for the downstream reaction. 

The velocity vector field and relative isobar distribution 
are shown in Fig.12 and Fig. 13 respectively. Because of the 
thermal expansion of the combustion products, the flame zone 
generates a pressure plateau region which retards and/or deflects 
the flow in front of the flame and accelerates the product gas 


behind the reaction front. Behind the flame tip another pressure 
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plateau occurring around the leading edge due to the presence of 
the stagnation point. The flow near the leading edge is strongly 
deflected by the combined contributions from the stagnation pres- 
sure and the shear stress induced by the vertical blowing velocity 
at the wall. Downstream of the leading edge, the magnitude of the 
pressure gradient is gradually decreasing and, at large x, a 
boundary layer type velocity profile is beginning to appear. 

The nondimensional blowing velocity (v,) distribution on the 
solid fuel plate is shown in Fig. 14. The blowing velocity is 
directly proportional to the solid fuel burning rate. The curves 
for Cases 1 and 2 are indistinguishable in this figure. For Case 
3, there is a finite, but still very small difference in blowing 
velocity compared with those of Cases 1 and 2. All these three 
cases are for envelope flames. The strongest blowing velocity in 
these cases occurs at the leading edge since the flame is closest 
to the fuel plate, and its magnitude is finite (v, = 0.644 in Cases 
1 and 2, v, = 0.612 in Case 3) instead of infinity as predicated by 
the boundary layer solution. The result given by the present 
theory is more realistic. As expected, the fuel evaporation rate 
drops quickly within a couple of thermal lengths, but continues to 
decrease with a more gradual rate in the downstream. More detailed 


comparison of the burning rate predictions between this study and 
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that using boundary layer fast kinetic model can be found in 


Section 4.6. 


4.3. Transition Between Envelope and Open-tip Flames 

As Damkohler number is decreased below certain value, the 
flame tip is found to become stabilized on the side of the fuel 
plate (open-tip flame). The draesyR ire from Peco thr i to open-tip 
flames can be apey by the flame tip isotherms in Fig. 15. 
Comparing Case 3 with Case 2, both envelope flames, although the 
Damkohler numbers are varied by a four-fold, the variation of 
thermal structure is insignificant. On the other hand, a slight 
decrease of Damkohler number from Case 3 to Case 4, we see that the 
flame front is pushed back at the leading edge region in Case 4, 
and the flame retreat process continues for Cases 5 sonnos In Case 
6, the flame tip becomes stabilized on the side of the fuel plate. 
Its location is more than two thermal distances downstream of the 
plate leading edge (x = y = 0). As a consequence, the fuel 
leading edge can no longer feel the heat feedback from the flame--a 
side-stabilized flame is established. 

The sequence of isotherm retreat in Fig. 15 can also be used 
to study the opening-up process for the envelope flame tip. For 
€xample, in Case 3, the isotherms for T= 5 and T= 6 reach the 
symmetry line (y = 0), but in Case 4 they do not. For lower 


Damkohler number, Case 5 shows this process is accelerated and now 
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the isotherm T = 3 does not reach y = o, but its apex is still 
ahead of plate leading edge. Finally, the apex of isotherm T = 3 
shifts behind the leading edge in Case 6 which is the end of 
transition. Clearly, we are observing a process where the envelope 
flame tip is broken up. 

The transition process between the envelope and the open-tip 
flames apparently LaveT¥es the shortening of the flame standoff 
distance, an increase of heat loss 0 the fuel plate and a 
temperature decrease in the neighborhood around the symmetry line 
(y = 0) in front of the leading edge. The transition is a 
continuous process, but occurs within a relative narrow range of 
Damkohler numbers. 

To define more precisely this range of transition, the fol- 
lowing criterion is were. The Hepinerie £987 the transition is 
marked when the surface blowing velocity at the leading edge (x = y 
=o) is at 95% of the maximum value (v, = 0.644) and the end of tke 
transition is taken when the leading-edge blowing velocity drops to 
zero. Using Fig. 14, the range of transitional Damkohler numbers 
has been determined to be within 5.82 x 10° and 7.03 x 105. The 


former value, corresponding to v_ = 0 at the leading edge, is found 


v 
by interpolation. 


The computed flame profiles for the two transitional flames 


(Cases 4 and 5) are presented in Figs. 33-42. 
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4.4. Open-Tip Flame 

Two cases of side-stabilized open-tip flames have been 
computed: Case 6 (Da = 5.56 x 105) and Case 7 (Da = 5.32 x 10°). 
The computed flame profiles for Case 6 are given in Fig. 16-20 to 
demonstrate the features of open-tip flames. 

Fig. 16 is the distribution of isotherms, A flame is seen to 
anchor itself on the sides of the fuel plate. The flame tip is 
located approximately at x = 2, which is also the location of 
maximum burning rate as shown in Fig. 14 (x = 2.2). For the side- 
anchored flame, there exists a low-temperature layer in the flame 
front region near the fuel plate (o <y <1) as shown in Fig. 16, 
which is resulted from the quenching effect of the cold fuel 
plate. The chemical reactivity is greatly reduced there and, as a 
consequence, both oxidizer and fuel vapor can leak through this 
quench zone. Fig. 17 is the distributions of mass fraction of fvel 
and oxidizer. It shows that the oxidizer penetration is especially 
pronounced as it is reinforced by the convective flow. In the 
neighborhood of the quench zone, the fuel and oxidizer is partially 
premixed due to the leakage. 

The nondimensional fuel reaction rate contours are shown in 
Fig. l18where the region of maximum reactivity islocated about 
two thermal lengths above the fuel surface. This zone serves as 
the stabilizer for the continuation of downstream reactions. A 


low: reactivity layer near the cold fuel surface is apparent. The 
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reactive zone in the flame front is broadened due to the extended 
overlap of the fuel and oxidizer there. For side-stabilized 
flames, the shapes of the reactivity contours are found to be 
similar to each other except their locations are different, as 
Damkohler numbers are decreased, the x-positions of the flame tip 
are shifted downstream and the y-locations for the maximum re- 
activity are slightly raised. This can be illustrated in Table 
4.4, which lists the locations of the most reactive point for side- 
stabilized flames and also includes that of envelope flames. 

Fijg2 + 108s? *the “velocrtty vector field and Fig. 20) iserue 
distribution of relative pressure contours. Since no flames exists 
ahead of the leading edge, there is only one pressure plateau 
region which is centered on the stagnation point. However, the flow 
downstream of pressure plateau is strongly influenced by the flane 
front. The incoming stream flows inward and then outward due to 
the flow retardation and then is pushed aside by the inverse 
pressure gradient originated from the leading edge. The downstream 
flow in front of the flame is first decelerated and/or 
deflected by the thermal expansion and is then accelerated behind 
the flame. 

For Case 7 (Da = 5.32 x 10°), the flame tip is stabilized at 
x = 8.20. Further decrease of Damkohler number pushes the flame 
downstream out of the leading edge region (0 ¢< x < 10). The stable 


locations of the open-tip flames are therefore very sensitive to 


+1 


TABLE 4.4 


Location of the Most Reactive Point 


Case Da Type of Flame Me Bas A 
1 4.50x107 E* upp wy 
2 Bi xi Oe E -2.40 0.795 
3 7.03x10> E RAO 0.795 
4 6.38x10° T -0.44 Oi Fe 
5 6.08x10> | T 0.64 2.09 
6 SES GR1 OF O ee DEO 
7 532x107 O BO 3.45 


* E: Enveloped Flame 
T: Transitional Flame 
O: Open-tip Flame 
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the magnitude of Damkohler numbers. 


4.5. Similarity Consideration for Envelope Flames 

Because of the use of nondimensional analysis, similarity 
relationships can easily be established. Exact similarity requires 
that all the nondimensional groups be the same, and this includes 
Damkohler number. However, in the present study, it is discovered 
that a number of flame characteristics are relatively insensitive 
to the variation of Damkohler number for envelope flames. In other 
words, for these flame properties, an approximate similitude exists 
independent of Da, as long as Damkohler numbers are large enough 
such that envelope flames are present. 

From computations for Cases 1-3, the following nondimensional 
quantities are found to be approximately invariant to Da: fuel 
burning rate (blowing velocity), locus of maximum temperature 
(flame envelope or flame standoff distance) and velocity field. In 
addition, temperature and species are also approximately the same 
except near the highest reactivity region. 

The small variation of the nondimensional blowing velocity 
from Case 1 to Case 3 has been shown in Fig. 14. Fig. 21 compares 
the locus of maximum temperatures. This locus is defined as the 
flame stand-off distance and is drawn from positions of maximun 
temperature along straight lines issued from the origin (x = y = 0) 


in the various directions. The largest discrepancy of the nondi- 
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mensional flame stand-off distance occurs at the stagnation stream- 
line which is about 17%. This covers over a range of Darkohler 
numbers from 4.5 x 10’ to 7.03 x 10°, so, indeed, an approximate 
independence from Da is observed. 

Damkohler number can be varied by changing either the fre- 
quency factor B or the free stream velocity U,, Eq. (2.19). 
Varying B without changing other parameters is difficult to achieve 
in practice, so we will focus on the influence of the free streazri 
velocity which enters into the nondimensional parameters only 
through Da. The Darkohler number is inversely proportional to U.?. 
Independence from the Damkohler number variation in the nondimen- 
sional results implies that the dimensional fuel burning rate is 
proportional to the free-stream velocity U, and the flame stand- 
off distance is inversely proportional to U,. In the limit of 
small convective velocity, a kinetically-strong (large Da), large- 
size but low-power (small burning rate) flame exists. This con- 
clusion is valid only in the absence of significant radiative heat 
losses [19]. 

Although a number of interesting quantities show approximate 
independence from the Damkohler number when they are expressed in 
nondimensional form, there are others which do not. The most 
conspicuous is the reaction rate. Figs. 11, 25 and 30 show that, 
not only the value of the maximum nondimensional reactivity 


decreases, but also the reactive domain broadens as TamkEohler 
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mumber is decreased. This is, of course, consistent with tke 


physical meaning of the Damkohler number. 


4.6. Relationship to Boundary Layer Combustion Model 

We mentioned in the introduction the necessity of elliptic-— 
type treatment in the leading flame front region when finite rate 
chemical kinetics are employed in the mathematical model. This 
region is critical for flame stabilization analysis which has been 
demonstrated in the results just presented. Downstream of the 
stabilization zone, boundary layer type flow is expected. Tke 
prediction from the stabilization zone analysis should provide the 
initial (upstream) conditions for downstream combustion calculation 
in the boundary layer. This section is devoted to studying the 
position and the process of transition from the elliptic to the 
parabolic regimes. 

The motivation to study the transition is a practical one. 
Computation for elliptic system (Navier-Stokes’ equations) is 
costly. When the computational domain or the number of grid points 
increases, the number of iterations needed for convergence rises 
quickly. Boundary layer type calculation, on the other hand, is 
fast by comparison since one space marching sweep is all that is 
required. The domain of elliptic region used in the computation, 
therefore, should be kept as small as possible. 


Since the transition from the elliptic region to the 


=} 


parabolic region is a gradual one, we list the profiles of U and T 
and their gradients at several x-positions (15, 20, 30). This is 
based on a computation for Case 2 parameters in which the boundary 


condition, Eq. (2.31), is applied at x,,, = 40. 


ma 
Table 4.5-1 shows that the u-velocity at x = 15 has a maximum 
at y = 12 and a local minimum (valley) at y = 25. The maximum 
Occurs just downstream of the heat release zone (flame) where the 
pressure field produced by thermal expansion accelerates the flow. 
The minimum occurs in front of the flame - a result of flow 
deceleration by the same mechanism. This peak-valley phenomena 
persists to downstream positions (x = 20 and 30) as can be seen in 
Tables 4.6-1 and 4.7-1, and is an indicator of the flame 
Oobliqueness. When the flame is exactly parallel to the flat 
plate, the velocity peak and valley should disapear which coincides 
with the boundary layer approximation. The velocity gradient du/dz 
is much smaller than unity at all x positions and the maximum 
value of the velocity gradient du/dy is of the order unity which 
indicate that the boundary condition, Eq. (2.31), taken at x,,, 
15 is a good approximation. The ratio of x- to y- velocity gra- 
dients is small close to the plate (y < 10) but has two infinite 
peaks at positions of vanishing du/dy (peak and valley). On the 


other hand, when compared with (du/dy) the ratio in the last 


Fines 


column is always much smaller than unity. In general, the approzxi- 


mation du/dx << @u/dy used in the boundary layer theory is con- 
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TABLE 4.5-1 


u(y) profile and its derivatives at x=15 
= ea Epspte leek lee 
y i ax dy 3x 7 sy) 0 
0.0000 0.0000 0.0000 0.3450 0.0000 0.0000 
0.1000 0.0345 0.0015 0.3214 0.0048 0.0045 
OFS1 00 0.1020 QO700s1 0.3074 0.0100 0.0089 
0.5410 Oe 730 0.0049 0.2893 0.0170 0.0142 
0.7951 0.2465 0.0071 O.a7a7 0.0258 0.0205 
1.0746 Ona 0.0086 0.2569 py jE how 0.0249 
i asses 0.4020 0.0104 a he 0.0456 O..0303 
164203 0.4830 Ebi 5 Uae 0.2244 0.0547 0.0356 
~ O823 0.5665 0.0138 0.2065 0.0669 0.0401 
Page wh 8 be! 0.6519 0.0154 0.1899 0.0808 0.0445 
ao ond) eee 0.735635 0.0169 0.1706 0.0990 0.0490 
3.4469 0.8210 CO. Oly s 0.1542 Osis 0.0499 
3.9915 0.7050 0.0178 O,iPSS2 On tit 8 0.0516 
4, 3707 0.9860 0.0178 eal ot 0 sila 0.0516 
5.2498 14.0625 0.0169 0.0993 0.1701 0.0490 
5.9747 1.1345 0.0163 0.0803 0.2028 0.0472 
6. tt LaiSes 0.0144 0.0610 .2oe7 0.0418 
7.6494 ieee Ja OF Oi 0.0446 0.2481 0.0320 
8.6144 PN be, 0.0080 OFO028s 8 Sarg el = Of0zst 
9.67358 EES la 0.0043 Ono L os Ofoas7 OFOl1 ao 
10.8434 1.3405 0.0003 0.0000 0.0000 0.0009 
12.1278 1.3405 -0,.0046 -0. 0209 0.2206 0.0134 
3.5405 lee a0 -0.0154 -0.0499 0.3079 0.0445 
15.0946 reas -~O, 000 -0.0773 0.4319 0.0970 
16.8040 1.1010 -0. 0559 -0.0745 On7 207 0.1620 
18.6844 0.9610 -0.0479 =0.0297 1.6114 0.1389 
20 ia 0.8995 -0..0316 -0.0143 2.2146 0.0917 
23-0282 0.8670 -0.0209 -0.0038 Souls 0.0605 
Primi ys LD) 0.8575 -0.0132 0.0009 14.5429 0.0383 
28.2841 0.8600 -0.0098 0.0041 2.2810 0.0285 
1 eid 0.8725 -0.0052 0.0072 0.7247 0.0151 
34.6438 O.o760 0.0034 0.0085 0.3993 0.0098 
38.3081 Os927a 0.0095 0.0074 P2792 Oe as 
42.3390 18 ie jhe ew 0.0107 0.0052 2.0722 0. 0312 
46.7729 0.9805 0.0089 0.0027 3.3415 © by 8 aban =| 
ain Coe Oe 79S0 0.0046 0.0012 3.0023 0.0134 
te Lok 1.0000 0.0000 0.0000 2 2 © © & | 0.0000 
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0.0000 
0.1000 
0.3100 
0.5410 
0. 7951 
1.0746 
}. S6ail 
1.7203 
BRAS 
ae ~OUS 
aa PRS 
3.4469 
A od? 
4.5907 
5.2498 
3.9747 
6.7722 
7.6494 
8.6144 
7.6/2 
10.8434 
les 1278 
13.5405 
15.0946 
16.8040 
18.6844 
£0,752 
23.0282 
23-310 
28.2841 
co Ero b prs 
34.6438 
38.3081 
42.3390 
46.7729 
31.6502 
2iaQioe 


u 


0.0000 
0.0285 
0.0845 
651435 
0.2060 
On27 70 
0.3390 
0.4095 
0.4830 
0.5590 
0.6360 
0.7150 
0.7945 


a= OS tan 


ee 
1.0270 
1.0990 
1.1665 

Py on 
t.,2/795 
1.3205 
1.3500 
1.3650 
1.3615 
1.3385 
1.2769 
1.1680 
1.0415 
0.9785 
0.9450 
0.9305 
0.9240 
0. 9220 
0.9260 
0.9390 
&. 9630 
O. 2976 


a 


TABLE 4.7-1 


gu 
ox 


0.0000 
0.0007 
0.0015 
0.0026 
0.0037 
0.0045 
0.0053 
0.0062 
0.0070 
0.0078 
0.0082 
0.0086 
0.0089 
0.0089 
0.0087 


0.0078 
0.0073 
0.0062 
0.0050 
0.0034 
0.0014 
=O. 0012 
-0.0054 
=O. OCI9 
-0.0191 
-0.0303 
-0.0327 
=0 0181 
-0.0111 
-0.0067 
-0.0043 
-0.0027 
-0.0012 
0.0000 
0.0010 
0.0018 


u(y) profile and its derivatives 


at x=30 
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TABLE 4.6-1 


u(y) profile and its derivatives at x=20 


gu ou poe) j 27 (22) 
y be 3x ay gx’ dy ax’ oy 
0.0000 0.0000 0.0000 0.3050 0.0000 0.0000 
0.1000 0.0305 0.0005 LPs ba Ws 0.0017 0. O@ie7 
0.3100 0.0925 -0.0019 ~2814 0.0069 0.0064 
0.5310 Oe laieie 0.0034 0.2656 0.0127 0 .O,2ie 
O07. 7951 0.2250 0.0045 D ehimleons 0.01/79 0.0148 
1.0746 ON are 0.0058 (oprenhche | 0.0244 0.0191 
1.3821 0.3690 0.0072 0.2247 0 Osis 0. Ona 
p eer AS oC} 0. 4450 0.0084 0.2097 0.0400 0 .Ogiao 
2.0923 8 wr { 8) 0.0094 0.1967 0.0479 0.0309 
Pad Bp. Ves} 0.6035 0.0105 0.1799 0.0586 0.0346 
a SSt7 0.6845 0.0116 0.1656 0.0699 0.0379 
3.4469 0.7665 OVOls2 6 A has fre 0.0809 0.0399 
3.9915 0.8485 0.0128 0.1344 0.0982 0.0420 
4.5907 0”. 9290 On Ose On LLL O. L1G0 0.0430 
$.2498 9; ls ae 0.0132 0.1028 Oe 0.0433 
5.9747 1.0820 O.012 0.0865 0.1491 0.0423 
6. 7/7 au ee eh ee 0.0123 0.0707 0.1738 0.0403 
7.6494 Da he «) O.Oilst 0.0554 0.1994 0.0362 
8.6144 1.2665 0.0095 0.0396 0.2406 0.032 
9.67358 1.3085 0.0073 0.0240 O's Oe 0.0238 
10.8434 1.3365 0.0042 Ae © lh vA 0.3594 0.0138 
ibs BPEL Mooole 0.0007 -0.0028 ks pedal 0.0023 
13.5405 1.3475 -0.0038 -0.0158 0.2403 0.0124 
15.0946 1.3230 -0.0094 -0.0433 0.2176 0.0309 
16.8040 1.2490 -0.0211 -0.0614 0.3433 0.0691 
18.6844 1 ee EX pe -0.0333 -0.0643 0.5174 0.1098 
20.7529 1.0005 -0.0406 -0.0284 1.5412 0.1333 
23.0282 0.9405 -0.0251 -0.0130 1.9315 0.0822 
dpi ae ke fp 0.9080 -0.0176 -0.0047 3.7290 0.0577 
28.2841 0.8950 -0.0123 -0.0005 24.8017 0.0403 
Bleacien 0.8935 -0.0089 0.0015 5.9338 0.0292 
34.6438 0.8985 -0.0062 0.0025 2.4087 0. O2.05 
38.3081 0.9080 -0.0037 0.0047 0.7818 . OL 24 
42.3390 0.9270 0.0002 0.00353 0.0386 Q.0007 
46.7729 0.9505 0. O02! 2 AN OW heer 0.4194 0.0070 
31.6502 OF Se 0.002 0.0045 0.6179 O. 009 5s 
wide lt hes 0.99938 0.0007 0.0001 8.4588 GO. O02 
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sidered reasonable at all these x positions, but the approximation 
for pressure (p = constant for flat plate) fails there because the 
flame is yet far from parallel to the plate. 

The temperature data given in Tables 4.5-2, 4.6-2, 4.7-2 give 
a consistent picture. The x-temperature gradient is in general 
small except just downstream of the flame where its absolute 
Magnitude may reach 0.5, The ratio of dT/@x to (8T/3y), =. oeararte 
last column also indicates a local breakdown of boundary layer 
assumption. It appears that the approximation dT/d@y >>aT/dx is 
more stringent than the corresponding one for u-velocity. 

Far downstream of the leading edge, not only boundary layer 
flow should prevail, but also similarity solution should be 
Obtained. Self-similarity requires infinite fast chemical kinetics 
which is asymptically reached since the nondimensional reaction 
rate in boundary layer coordinates increases with downstream 
distance from the leading edge. Fig. 22 gives a comparison of the 
fuel blowing velocity in the near-leading edge region between the 
self-similar boundary layer solution and the present full Navier- 
Stokes’ solution. The boundary layer solution is obtained from 


Ref. 3 where the mass transfer number, By is defined by 


Be (Cs (Tosa ta) 0+ Os (rye ero) VL 


(C,(T. - Ty) + Q,Yoa]/L 


1.2896 
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0.0000 
0.1000 
0.3100 
0.5410 
Oa7951 
1.0746 
1 23621 
lie 7203 
2.0923 
2.9015 
2a7at7 
3.4469 
3 otter ti 
4.5907 
37-2498 
Ser nas 
6.77/22 
7.5494 
8.6144 
9.6728 

10.8434 

12.1275 

13.5405 

15.0946 

16.8040 

18.6844 

20. 4a? 

23.0282 

2o.aa10 

28.2841 

3l.312e 

34.6438 

38.3081 

42.3390 

46.7729 

Slee wOd 

7 6 Oli 


oe Cad] 


ay 


2.1665 
2.2490 
2.4180 
2.6000 
2.7910 
me 7200 
3.2020 
3.4280 
3.6660 
3.9130 
4.1690 
4.4370 
4.7100 
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TABLE 4.5-2 


ar 
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0.0025 
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TABLE 4.6-2 
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3.5660 
3.7830 
4.0130 
4.2510 
4.4980 
4.73540 
5- 0180 
a= ee) 
JeUwlO 
53-8240 
6.0930 
6.3250 
6 i 7 
6.8160 
7.90110 
6.9980 
5.3340 
3.3280 
1.8160 
1.3130 
1.1050 
1.0310 
1.0100 
1.0010 
1.0010 
1.0010 
1.0010 
1.0010 
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TABLE 4.7-2 


at 
ay 


0.6500 
0.6190 
0.5974 
0.5982 
0.35724 
0.5659 
0.5381 
0.5242 
0.4961 
0.4821 
0.4645 
0.4369 
0.4122 
0.5384 
0.35641 
0.3373 
0.35009 
0 aoa 
0.2534 


‘gp Js Le! 
er 


0.1853 
0.1593 
Os12S5 
-0.0076 
-0.8849 
-0.9698 
-0.6645 
-0.2010 
-0.0756 
-0.0244 
-0.0063 
=9 0025 
0.0000 
0.0000 
0.0000 
0.0000 
0.0000 


profile and its derivatives at x=39 


ar ,at 
beter 


0.0002 
0.0030 
0.0075 
0.90133 
0.0200 
0.0264 
0.0346 
0.0426 


‘Qe £0537 


0.0607 
0.0705 
0.0834 
0.0968 
0.1096 
0.1220 
0.1396 
0.1626 
0.11759 
OS1ISS 
0.2354 
0.2828 
0.3175 
0.3818 
wei OFT, 
0.0009 
0.3642 
0.6421 
1.5801 
1.3984 
1.7093 
2.4034 
1.7506 
KREXKE 
KXXEXX 
-2 2 2S ® | 
XXKXXX 
KXKKXE 
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0.0028 
0.0069 
0/Otas 
0.0176 
0.90230 
0.0287 
0.0343 
0.0409 
0.0450 
0.0504 
0.0561 
0.0614 
0.0655 
0.0684 
0.0724 
0.0753 
0.0765 
0/0/78 
0.0806 
0.0806 
0.0778 
00737 
0.0668 
0.0013 
0.5434 
0.6564 
0.4885 
0.162 

0.0643 
0.0233 
0.0066 
0.0028 
0.0000 
0.0000 
0.0000 
0.0000 
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Using this value for B., the blowing velocity for the boundary 


layer is found to be 
vy = 0.1143/x1/2 


This is plotted in Fig. 22 to compare with the result of the 
present calculation. Exceptat the very leading edge region (zx 
0.05), the blowing velocity from the Navier-Stokes solution is 
higher. This is consistent with past findings on skin friction 
calculation and measurement [27, 28]. The present computation 
suggests that it may take a long distance (in terms of thermal 
length) in order for the blowing velocity to drop and reach the 
asympotic value predicated by the self-similar solution. This 
should not te too surprising since tke downstream of the leading 
edge the boundary layer thickness is many times larger than the 
thermal length. 

While we tkink for the purpose of analyzing flame 
stabilization, it is adequate to use a reasonable small 
Computational domain, it is not clear at this stage the degree of 
inaccuracy induced if we apply the initial conditions for the 
boundary layer at the downstream location of the elliptic region 
used in the present work. In principle, boundary layer requires 
large local Reynold number. In this study the highest Reynold 
number is only 40. More detailed work in this respect appears to 


be needed. 
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4.7. Comparisons With Experiments and Other Related Studies 

To our knowledge, there is no experiment corresponding ex- 
actly to the model configuration analyzed here. The comparison 
with experiments will have to be on limited aspects of results, and 
is qualitative and fragmented in nature. This applies also to the 
comparison with other theoretical works. 

Configuration-wise, the experiment closest to the theory 
appears to be Chu’s work [20] where an upward propagating flame 
Over a thin paper sample is investigated. The width of the paper 
sample he used, however, was very narrow and three-dimensional 
effect was obvious. In addition, the flow is naturally-induced in 
the experiment rather than forced such as assumed in the model 
Nevertheless, the temperature profiles measured at the center line 
bears much resemblance to the comptted profiles for the 
transitional flames (Cases 4 and 5). 

The computed temperature profiles for the open-tip flames 
agree qualitatively with measurements in Ref. 21 where a flame 
Spreads over a paper sample in an opposed natural convective flow. 
In particular, the measured temperatures downstream of the 
Pyrolysis front show fuller profiles on the fuel-rich side 
consistent with the present calculated results as well as with 
those in Refs. 10, 12 and 22. 

The overlapping feature of fuel and oxidizer near the flame 


tip quench zone for side-stabilized flames agrees with predictions 
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in Refs. 10 and 12. The computed reactivity contours are similar 
to those calculated in Refs. 12 and 18, 

The retardation and deflection of flow in front of the flame 
and flow acceleration behind the flame has been observed experimen- 
tally by Hirano et al. [23]. All of the Navier-Stokes’ calcula- 
tions allowing gaseous thermal expansion [Refs. il and 12 and 
present work] reproduce this phenomena. The retreat of the flame 
tip for side-stabilized flames has also been found in Ref. 11 for 
premixed flames and in Ref. 13 for diffusion flames. Fxrperimen- 
tally, the observation by Ray [24] is also consistent with the 
model finding. Ray observed that in opposed-flow flame spreading 
experiment in a wind tunnel [25], a flame could become quenched as 
it propagated toward upstream. Although the mean velocity in the 
tunnel is kept constant, the flame sees a thinner boundary layer as 
it spreads upstream. Thus, a velocity profile effect is observed. 

We are not aware of any experiments on envelope flame blowoff 
at the leading edge of the fuel plate. The classical experiment by 
Spalding [26] using a porous sphere did show that as the forced 
velocity was increased, the envelope flame broke up and became 
anchored on the sides. This similarity with our calculation is 
€ncouraging, but it has to be recognized that the two geometries 
are quite different. In general, there is a lack of experimental 


data on multi-dimensional envelope diffusion flames. 


CHAPTER V 
CONCLUSION 


A model of a diffusion flame at the leading edge of a fuel 
plate in a forced convective flow has been developed and solved 
numerically to study the flame stabilization and blowoff phenomena. 
In addition to finite-rate chemical kinetics, the streamwise heat 
Conduction and mass diffusion are included as the essential 
mechanisms for the reaction initiation in the flame tip 
Stabilization zone. The mathematical system consists of two 
dimensional Navier-Stokes’ momentum, energy and species equations. 
A one-step overall chemical reaction with a second-order, finite- 
rate Arrhenius kinetics is adopted. 

A nondimensionalization procedure is performed to identify 
the dimensionless parameters, ten have been found. The numerical 
Tesults presented, however, are restricted to the effects of 
varying Damkohler number, which is proportional to the ratio of 
flow residence time in one thermal length to chemical reaction 
time. 

When Damkohler number is large, envelope flames surrounding 
the plate leading edge are found. As Damkohler number is reduced, 


flame retreats and an open-tip flame stabilized on the side of the 
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fuel plate can be achieved. The position of the flame tip is very 
sensitive to the magnitude of the Damkohler number. Transition 
between envelope and open-tip flames occurs in a narrow range of 
Damkohler numbers. The flame structure for all the three stages 
are presented graphically. 

Contrary to side-stabilized flames, envelope flames have been 
found to be insensitive to Damkohler number for many of their flame 
characteristics. This observation has been exploited, especially 
in relation to the effect of varying the free stream velocity. In 
the absence of radiative heat losses, a kinetically-strong (large 
Da), large-size, but low-power (small burning rate) flame exists in 
the limit of small convective velocity. 

The use of the present model to provide the upstream 
(initial) conditions for boundary layer combustion calculation was 
discussed. There are many other applications of the model, 
Particularly in the area of diffusion flame spread and spread 


limits. These will be explored in the future. 


10% 
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APPENDIX A 


DERIVATION OF FINITE DIFFERENCE EQUATIONS 


A.1 Introduction 
The momentum, energy and species equations (Eqs. (2.20) - 
(2.26)), with the continuity equation, can be rewritten in the 


following general form: 


3 36 nae _ peo) 2g (A.1) 
a chase (WAVES, © SY) .Tar9 
where @ is the dependent variable, [ and S 4 represent chewcon— 


ductance and the source terms relevant to >. Their expressions 
are listed on Table A.1l. 

The finite difference equation for each dependent variable is 
derived by integrating Eq. (A.1) over the corresponding 
computational cell introduced in Chapter III, the integrations of 
temperature, fuel and oxidizer mass fractions are performed over 
the main cell, the integrations for the velocity components u and 


v, are carried out over the u cell and the v cell respectively. 
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The resulting terms are expressed through the grid-point values. 


The integrating procedures for all these equations are the 
same. Therefore, a detailed derivation for the variable 
integrating over the main cell is sufficient to illustrate the 
procedure and it is given in the next section. A summation of the 
results of integration for all of the variables is presented in 


Section A.3. 


A.2. Integration Procedure 
In the main cell, shown in Fig. 3, attention is focused on 
the grid point P, which has the points E, W, N and S as its 
ares The letters e, w, n and s denote the faces of control- 
volume which are located at the mid-way between the grid points. 
The integration of Eq. (A.1) over the main cell (control- 


volume) is 


(A.2) 
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The first integral on the left-hand side of Eq. (A.2) can be 


written as 


n <5 n 30,¢e 
Is sx (pud - ree 2) dxdy = C [oud - Toyd, dy 

Ss 

= I, mie ae 
(A.3) 
where 
n 
rE . rad 

Ir, = [. [pug xe dy, 

I, = li [pud - rao), rob tae 

Approximately, I, is given by 

ad LA MeL) 2 tcp Ah 2) le 

e a axe ew (A.4) 


ue 


where See = (bye5 ~ Ayso) is the east or the west face area 
of the control volume. Eq. (A.4) consists of convection and con- 
duction/diffusion terms. The most straightforward way to obtain 
the expression for both terms is assuming a piecewise-linear pro- 
file (central difference scheme) between the grid points. However, 
this kind of treatment sometimes causes disastrous outcome such 
that the finite difference equation may diverge. Therefore, in 
this work upwind scheme is adopted which has been first put forward 
by Courant, Isaacson and Rees (1952). The upwind scheme is applied 


on the convection term only, whereas the conduction term is still 


using the central difference scheme. Using this, Eq. (A.4) becomes 


ed Eopec rete Dobe ne Ea. (A.5a) 
T, = [2(O, + Pp)up% 2 eae 
PF 
rs up? 0, 
or 
ree ad oe (0) 
1 Pp oases Pi. 
I. = C2(P, + Pp) Updos =| SoMa Exe, ew (A.5b) 
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Let 
ty) 
me = £(0p + Op)UDAL, 
nari Bip 
o P E,,9% 11-0 
ds = aaPEY JAcw * £|m3 |. 
PE 


(A.6) 


With the definitions given by Eq. (A.6), the conditional statements 


(A.5) can be written in a more compact form: 


> 1-9 > 1.9 
I. = (d, + zm 14, - [d. - zmp]on. (A.7) 


Using a similar manipulation for I, in Eq. (A.3), we have 


(A.8) 


a 


where 
o ; a 
ene, a z(Py . Pp) Uyraw! 
Le eae 
ts) W P,,% 7) 
dq. = 4( 7 = HF as * z|m, é 
*WP 


Substituting Eqs. (A.7) and (A.8) into Eq. (A.3) we get 


ne 
emer et 
| | Fx (pud - [a5 ) dxdy 


s‘w 
= Tr, ~ I, (A.9) 
) ) ) ) d a) E > 1 
Similarly, the second integral on the left-hand side of Eq. 
(A.2) can be written as 
e-rn 
0 _ roe 
| | aplove Pay dydx (A.10) 
_ o 7) ) 1_ 9 ) 1_ > d , > 
= Cd. + 3m, + d< - 2m.1%, = Cd. - mm 19, ~ [dq +2mo]o, 
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where 
> 1 d 
Sra ) 
Ms = 2(95 + P_)VCAL 
Ee Sod 
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Dean eS Ie Bk 11% 
seh (Atparycror ace zimol, 
Ysp 
and a? = (Ax? + ax® ) is the north or south face area of 
ns ~ 2'°°PE WP 


the control volume (main cell). 

The resulting appearance of integration for the source term 
on the right-hand side of Eq. (A.2) depends on the actual 
expression Sy itself. However, the integrating procedure is the 
same as what has been presented. In order to achieve a converged 
solution, the resulting finite difference representation of the 


source term is suggested to be in the following form: 


(A.11) 


a9 


where the quantity : must not be positive to make numerical schene 
stable. 

Finally, substituting Eqs. (A.9), (A.10), and (A.11) into Eq. 
(A.2), the differential equation can be transformed into a standard 


form of fully implicit finite difference equation: 


ty) fo) b oO. 4 (i) 
>’ p = an?s “ anew + axon * a5%> * Sas (A.12) 
where 
«) 
a - de - 2M5, 
ae = a? + ime, 
an = ay = 2M 
ae - ae + ime, 
do  .d ) fs) ) ty) ) 
ap = ae + ay + an + ac + (m5 - m,+m, - ma) 
= Sp 
sd? . s? 
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The coefficients, ap, a,, a, and a,, are guaranteed to be positive 
quantities by the definition themselves. For convergence, the 
value of a° should be positive. A special treatment for the mass 
3 g. 789 od; Ay 
source term, (mg - my + my - 3g), inside a, is needed. The finite 


difference equation fox continuity equation is obtained by integra- 


ting Eq. (2.20) over the main cell. The result is 


e e ote E (A.13) 
Coup - Cyty + CnYp sys 
where 
Cop = $(0,%+ 70 ya? 
E =ap E’*ew! 
Cer lea (Dine nek 
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Ce she ((DLG4op ya® 
N P NaanSs. 
Cree i) 


If mass balance is achieved in the control volume, from Eq. 
(A.13), the term of (mp - my + my - mg) vanishes in Eq. (A.12). 


During the iteration procedure, however, this term usually is not 
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zero and may lead to computational instability. The trick used to 


avoid this condition is to replace ao and s? in Eq. (A.12) by 


d db > ) ) ) 0) d ) 
ap = ap + ay + ay * ag + Max(m, - Mm + My - Mor Oj] - Sp, 
(A.14) 
se = 34 + Max(me - m + my - me, 0105 (Al 1S) 


where superscript + denotes the quantity resulted from last 
iteration. This will ensure that all of the coefficients in Eq. 
(A.12) are positive to avoid the divergence due to the existence of 
negative coefficient in a0. 
A.3. Summation of Finite Difference Equations 

The forms of finite difference equations for T, Yr, and Y, 
are identical except the expression of some of the coefficients. 
The form of the equation is given by Eq. (A.12) except the 
coefficients ad and s¢ are replaced using Eqns. (A.14) and (A.15). 


The undefined [, 5? and Sy for the variables T, Yr and an are 


tabulated in Table A.2. 
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for 6= v, the finite difference equation is 
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APPENDIX B 
CONSIDERATION OF NUMERICAL ERRORS 


The accuracy of computational solutions depends on a number 
of parameters used in the numerical scheme, such as the grid mesh 
size, the finite location of boundary conditions and ne 
convergence criterion in the iteration procedure. Among these 
factors, the grid mesh size is of primary importance. Generally 
speaking, decreasing the mesh size increases the accuracy. But, 
the use of finer mesh size, which also increases the mesh number, 
requires more computational time and memory storage and is 
sometimes impractical. Compromising mesh size normally is 
determinec through numerical experiments. Tests are made with two 
sets of mesh points, where ranges are 49 x 38 and 69 x 48 
respectively, in the same computational domain described in Section 
3.6 for Case 2 (Da = 2.81 x 10%). In the case of mesh points 69 x 
48, we change the exponents from 1.10 to 1.05 in both +z and -x 
directions and to 1.068 in +y direction. Therefore, grid sizes 
for 69 x 48 are decreased averagely 40% in both +x and -x direc- 
tions and 25% in vidizece ion with respect to the 49 x 38 mesh. The 
comparison is plotted in Fig. Bl, which is the temperature dis- 


tribution T versus the vertical distance y in three streamwise 
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locations, x = -2.40 (upstream), x = 0 (leading edge) and x = 10 
(downstream). The largest difference occurs at x = 0, because the 
smallest grid size concentrates on that region. For this 
comparison, the mesh-point number 49 x 38 is acceptable since the 
profile in the fuel side of the flame is indistinguishable from the 
one with 69 x 48 mesh and the local relative shift of the 
temperature profile in the oxidizer side of the flame is less than 
10%, : 

An iterative process is said to have converged when further 
iteration will not produce any change in the values of the 
dependent variables. In practice, the iterative process is 
terminated when a convergence criterion is satisfied. A meaningful 
method of monitoring convergence is to examine how perfectly the 
finite difference equations are satisfied by the current values of 
the dependent variables. For each grid point, a residual R can be 


calculated from 


R= aap apy * So - an?» 

Obviously, when the discretization equation is completely 
satisfied, R will be zero. A suitable convergence criterion is to 
require the largest value of |R | be less than a certain small 


number, e. In this work, the quantity of e is set to be 0.005. 
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When the conditions of 


[Rr | max Se, amas Se, |R,| max < & 


are satisfied simultaneously for temperature T and both velocity 
components u and v, then the iteration procedure is stopped and the 
results are printed out. 

The number of iterations to obtain the one rents for 
envelope flame is about 120. Each iteration includes four 
alternative sweeps (ADI) in x and y directions for every dependent 
variable, and the corresponding CPU time spent is approximately 20 
minutes. For the side-stabilized flame, the iteration number is 
about 300, which is approximately one hour CPU time. The 


computational work was performed using a DEC VAX 11/780. 


APPENDIX C 


COMPUTER PROGRAMS 


(I) INTRODUCTION OF PROGRAMS 


ARE: 


THE COMPUTER CORE CONSISTS OF 13 PROGRAMS,AND THEY 


C1) 


ee 


(7) 


(8) 


MCELLCFORS 


UCELL. FOR: 


VCELL.FOR: 


TCELL.FOR: 


YCEED.EOR: 


MAIN FROGRAM TO SOLVE THE FRESSURE 
CORRECTION EQUATION, UPDATE THE 
BLOWING VELOCITY ON THE. FUEL PUaATe 
AND DECIDE WHEN THE PGOGRAM SHOULD 
STOF AND FRINTOUT THE RESULTS. 


SUBFROGRAM TO SOLVE X-MOMENTUM 
EQUATION. 


SUBFROGRAM TO SOLVE Y-MOMENTUM 
EQUATION. 


SUBPROGRAM TO SOLVE ENERGY EQU. 


SUBPROGRAM TO SOLVE SPECIES EQU: 


F1.FOR: SUBPROGRAM TO EVALUATE THE MASS FLUX 
AND DIFFUSION COEFFICIENTS AT THE 
CONTROL-VOLUME SURFACE IN X-DIRECTION 
FOR UCELL AND VCELL. 


F2.FOR: SUBPROGRAM TO EVALUATE THE MASS FLUX 
AND DIFFUSION COEFFICIENTS AT THE 
CONTROL-VOLUME SURFACE IN Y-DIRECTION 
FOR UCELL AND VCELL. 


F3.FOR: SUBF'ROGRAM TO EVALUATE THE MALL FLUX 
AND DIFF. COEF. AT THE CONTROL-VOLUME 
SURFACE IN BOTH X AND Y DIRECTION FOR 
TCELL AND YCELL. 


90 


(9) 


(10) 


Ot 


F4.FOR: SUBPROGRAM TO EVALUATE THE COEFS. OF 
ENERGY AND SPECIES DIFFERENCE EQS. 


SPU.FOR: SUBPROGRAM TO EVALUATE THE SOURCE 
TERN IN X-MOMENTUM EQUATION. 


SPV.FOR: SUBPROGRAM TO EVALUATE THE SOURCE 
TERM IN Y-MOMENTUM EQUATION. 


GRIDO.FOR: SUBFROGRAM FOR GRID GENERATION. 


TRAD.FOR: SUBPROGRAM OF TDMA. 


(II) INPUT DATO FILES 


C2) 


DI.DAT: CONTROLLING DATA FILE 
2222922222992 2292 2992399229289 9 0283 2: 
MI MJ MH MRUN MALL MWRPF MUV MTFO & 
RE PR LE TW TA YFA YOA x 
Datieech i tUCCRITV CRITT 
WI WJ WE WL WM 


x 
x 
WV WP K 
x 
x 


ot pe HS MH HS MH 


DA GF FO ET HLT 

2 2929229222 2922992229292 5323292 23 8 

WHERE 

MI: NO. OF GRID IN X DIRECTION 

MJ: NO. OF GRID IN Y DIRECTION 

MH: POSITION OF LEADING EDGE AT Y=0. 

MRUN: NO. OF ITERATION 

MALL: =1, ONLY READ INPUT HOT FPROFILES; 
DT.DAT,YF.DAT,YO.DAT. 
=OTHERS, READ ALL INFUT PROFILES, 
INCLUDING U,V,P,T>,YF>,YO. 

MWRF: =1, PRINTOUT REACTIVITY DATA 
=OTHERS, NO PRINTOUT FOR REACT. 

RE: VALUE OF REYNOLDS NO. 

PR: VALUE OF PRANDTL NO. 

LE: VALUE OF LEWIS NO. 

TW: FUEL-PLATE TEMPERATURE 

TA: AMBIENT TEMPERATURE 

YFA: AMBIENT FUEL MASS FRACTION 

YOA: AMETENT OXIDIZER MASS FRACTION 

DAMF: VALUE OF INERTIA FOR PRESSURE 
CORRECTION EQUATION. 

CRITU: CONVERGENCE CRITERION FOR U 

CRITV: CONVERGENCE CRITERION FOR V 

CRITT: CONVERGENCE CRITERION FOR T 
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WI: UNDER-RELAXATION FACTOR FOR U 

WJ: UNDER-RELAXATION FACTOR FOR V 

WK: UNDER-RELAXATION FACTOR FOR YF,YO 

WL: UNDER-RELAXATION FACTOR FOR T 

WM: UNDER-RELAXATION FACTOR FOR VWALL 

WV: UNDIER-RELAXATION FACTOR FOR VELO- 
CITY CORRECTION. 

WP: UNDER-RELAXATION FACTOR FOR PRE- 
SSURE CORRECTION. 

DA: VALUE OF DAMKOHLER NO. 

QP: HEAT OF REACTION PER UNIT MASS 
FUEL: 

FO: STOICHIOMETRIC MASS RATIO. 

ET: ACTIVATION ENERGY. 

HLT: LATENT HEAT. 


{2)5 HOT) PROFICES: 
CUT.DAT: INITIAL TEMPERATURE PROFILE 
YF.DAT: INITIAL FUEL MASS FRACTION PROFILE 
YO.DAT: INITIAL OXI. MASS FRACTION PROFILE 


(III) SEQUENCE OF EXECUTION 


(1) FOR MCELL 
FOR UCELL 
FOR VCELL 
FOR TCELL 
FOR YCELL 
FOR F1 
FOR F2 
FOR F3 
FOR F4 
FOR SFU 
FOR SPV 
FOR GRIL 
FOR TRAD 


(2) LINK MCELL,UCELL,VCELL, TCELL, YCELL, 
F1,F2,F3,F4,SPU,SPV,GRID, TRAD 


(3) ASS DI.DAT FOROO! 
ASS DT.DAT FOROOS 
ASS YF.DAT FOROO7 
ASS YO.DAT FOROOS 


(4) RUN 


OUTU. DAT 
OUTV. DAT 
OUTF. DAT 
OUTT. DAT 
OUTF. DAT 
OUTO. DAT 
OUTW. DAT 


MCELL 
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FORO10 
FORO11 
FORO12 
FORO13 
FORO14 
FORO1S 
FORO16 
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(IV) PROGRAMS! 
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Figure 3. Main cell 
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